function [score,dPyy,dInov,gradMom,errorMes] = getRowOfScore_mp(paramsValues,setupFilter,outFilter,epsValue,stepOption,pos_i)

[samples,dimy] = size(setupFilter.data);
score = zeros(samples-setupFilter.numInitial,1);
dPyy  = zeros(dimy,dimy,samples);
dInov = zeros(dimy,samples);

% Positive step
paramsValues_epsp = paramsValues;
paramsValues_epsp(pos_i) = paramsValues_epsp(pos_i) + epsValue;
[~,errorMes_epsp,~,~,outFilter_epsp] = objectFuncQML(paramsValues_epsp,setupFilter);
if errorMes_epsp == 1
    errorMes = 1;
    return
end

% Negative step
if stepOption == 2
    paramsValues_epsm = paramsValues;
    paramsValues_epsm(pos_i) = paramsValues_epsm(pos_i) - epsValue;
    [~,errorMes_epsm,~,~,outFilter_epsm] = objectFuncQML(paramsValues_epsm,setupFilter);
    if errorMes_epsm == 1
        errorMes = 1;
        return
    end
end

if stepOption == 1
    score(1:samples,1) = [zeros(setupFilter.numInitial,1);...
                         (outFilter_epsp.LogL(setupFilter.numInitial+1:end) ...
                        - outFilter.LogL(setupFilter.numInitial+1:end))/epsValue];
    for t=1:samples
        numObs = sum(~isnan(setupFilter.data(t,:)));
        dPyy(1:numObs,1:numObs,t,1) = (outFilter_epsp.Pyy(1:numObs,1:numObs,t) - outFilter.Pyy(1:numObs,1:numObs,t))/epsValue;
        dInov(1:numObs,t,1)         = (outFilter_epsp.inov(1:numObs,t)         - outFilter.inov(1:numObs,t))/epsValue;
    end
    if isfield(outFilter,'momCondition') == 1
        % The score of the moment conditions
        gradMom(:,1) = (outFilter_epsp.momCondition-outFilter.momCondition)/epsValue;
    end
elseif stepOption == 2
    score(setupFilter.numInitial+1:samples,1) = (outFilter_epsp.LogL(setupFilter.numInitial+1:end) - outFilter_epsm.LogL(setupFilter.numInitial+1:end))/(2*epsValue);
    for t=1:samples
        numObs = sum(~isnan(setupFilter.data(t,:)));
        dPyy(1:numObs,1:numObs,t,1) = (outFilter_epsp.Pyy(1:numObs,1:numObs,t) - outFilter_epsm.Pyy(1:numObs,1:numObs,t))/(2*epsValue);
        dInov(1:numObs,t,1)         = (outFilter_epsp.inov(1:numObs,t)         - outFilter_epsm.inov(1:numObs,t))/(2*epsValue);
    end
    if isfield(outFilter,'momCondition') == 1
        % The score of the moment conditions
        gradMom(:,1) = (outFilter_epsp.momCondition-outFilter_epsm.momCondition)/(2*epsValue); 
    end
end

end